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Abstract 

In this paper, we outline a model of graph (or network) dynamics based on two ingredients. 
The first ingredient is a Markov chain on the space of possible graphs. The second ingredient is a 
semi-Markov counting process of renewal type. The model consists in subordinating the Markov 
chain to the semi-Markov counting process. In simple words, this means that the chain transitions 
occur at random time instants called epochs. The model is quite rich and its possible connections 
with algebraic geometry are briefly discussed. Moreover, for the sake of simplicity, we focus on the 
space of undirected graphs with a fixed number of nodes. However, in an example, we present an 
interbank market model where it is meaningful to use directed graphs or even weighted graphs. 
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I. INTRODUCTION 



The publication of Collective dynamics of 'small world' networks by Watts and Strogatz 
Q gave origin to a plethora of papers on network structure and dynamics. The history of 
this scientific fashion is well summarized by Rick Durrett |2J: 

The theory of random graphs began in the late 1950s in several papers by Erdos 
and Renyi. In the late twentieth century, the notion of six degrees of separation, 
meaning that any two people on the planet can be connected by a short chain 
of people who know each other, inspired Strogatz and Watts [l| to define the 
small world random graph in which each side is connected to k close neighbors, 
but also has long-range connections. At about the same time, it was observed 
in human social and sexual networks and on the Internet that the number of 
of neighbors of an individual or computer has a power law distribution. This 



inspired Barabasi and Albert {3] to define the preferential attachment model, 
which has this properties. These two papers have led to an explosion of research. 
While this literature is extensive, many of the papers are based on simulations 
and nonrigorous arguments. 

Incidentally, the results of Watts and Strogatz were inspired by the empirical and theoret- 



ical work by Milgram (4] and Granovetter [5| back in the 1960s and 1970s; similarly, the 
preferential attachment model by Barabasi and Albert is closely related to the famous 1925 
paper by Yule js] as well as to a celebrated work by Herbert Simon published in 1955 [3] 
(see also chapters 8 and 9 in reference {§]] for a recent analysis on Simon's results). This 
body of literature is partially reviewed in Durrett's book [2] as well as in a popular science 
book written by Barabasi [9]. 

It might be interesting to understand why this scientific fashion was born and how. On 
this respect, we can quote Wikipedia's article (as retrieved on 4 May 2011) on Milgram's 
experiment in popular culture jlo| : 

Social networks pervade popular culture in the United States and elsewhere. In 
particular, the notion of six degrees has become part of the collective conscious- 
ness. Social networking websites such as Friendster, MySpace, XING, Orkut, 
Cyworld, Bebo, Facebook, and others have greatly increased the connectivity 



2 



of the online space through the application of social networking concepts. The 
"Six Degrees" Facebook application calculates the number of steps between any 
two members. [. . .] 

In other words, the social character of human beings combined with the hyper-simplification 
(trivialization) of some results promoted by leading science journals might have triggered 
interest in social networkology also outside scientific circles. Moreover, the emergence of 
social networks in the Internet has indeed made some tools developed by networkologists 
profitable. However, a deeper analysis by sociologists and historians of science will be 
necessary to falsify or corroborate such hypotheses. 

In this paper, we pay our tribute to this fashion, but we slightly depart from the bulk 
of literature on social network dynamics. First of all we consider time evolution also in 
continuous time and not only in discrete time. As the reader will see, this will be enough to 
give rise to interesting non stationarities as well as to non-trivial ergodic behavior. Moreover, 
to begin with a simple situation, we will be concerned with undirected graphs whose number 
of nodes M does not change in time. These restrictions can be easily overcome and, indeed, 
in the following, an example with directed graphs will be presented. The dynamic variable 
will be the topology of the graph. This approach is motivated by the following considerations. 
Social networks are intrinsically volatile. You can be in contact with someone for a finite 
time (at a meeting, during a phone call, etc.), but never meet this person again in the 
future. This interaction may or may not have effects on your future actions. If memory 
is not a major issue, the new configuration of the graph will only depend on the previous 
configuration. Memory is indeed an issue, but again, to simplify the analysis, we will consider 
a semi-Markov dynamics on the state space of all the possible graphs with M nodes. It is 
already quite rich. Incidentally, notice that, except for the case of infinite memory, finite 
memory processes in discrete time are Markov chains. 

The dynamics will be defined by a Markov chain subordinated to a generic counting 
process. Similar models have been around for many years. They were (and are) commonly 
used in engineering and decision analysis and, on this point, the interested reader can consult 
the monograph by Howard 

In this framework, it is often assumed that the waiting times between consecutive events 
do follow the exponential distribution, so that the corresponding counting process is Pois- 
son. Indeed, many counting processes with non-stationary and non-independent increments 
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converge to the Poisson process after a transient. If these counting processes are renewal, 
i.e. inter-arrival times { Ji}^i are independent and identically distributed (iid) random vari- 
ables, it is sufficient to assume that the expected value of these inter-arrival times is finite. 
However, recently, it has been shown that heavy-tailed distributed interarrival times (for 
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14] . After defining the 



which E(Jj) = oo) play an important role in human dynamics 
process in Section II, we will present two introductory examples in Section III and a detailed 
model of interbank market in Section IV. 

II. THEORY 

This section begins with the definition of the two basic ingredients of our model, namely 

1. a discrete-time Markov chain on the finite set of 2 M ( M+1 )/ 2 undirected graphs with M 
vertices (nodes), and 

2. a counting process N(t) for the point process corresponding to a renewal process. 
The rest of the section is devoted to the definition of the basic model class. 



A. Ingredient 1: a Markov chain on graphs 



Consider an undirected graph Gm — (Ym^E) where Vm represents a set of M vertices 
(nodes) and E the corresponding set of edges. Any such undirected graph can be represented 
by a symmetric M x M adjacency matrix Ag M , or simply A, with entries Aij = Aj^ = 1 
if vertices i and j are connected by an edge and A^j = Ajj = otherwise. Note that 
algebraic graph theory using linear algebra leads to many interesting results rela ting spe ctral 
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lij. For 



properties of adjacency matrixes to the properties of the corresponding graphs 
instance, the matrix 

A= l 1 
1 1 

corresponds to a graph where there are no self-connections and each vertex is connected to 
the other two vertices. As mentioned above, for a given value of M there are 2 M ( M+1 ^ 2 
possible graphs. To see that, it is sufficient to observe that the M diagonal entries can 
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assume either value 1 or value and the same is true for the M(M — l)/2 upper diagonal 
entries. Now, denote by Gm the set of 2 M ( M+1 )/ 2 undirected graphs with M nodes. Consider 
a sequence of random variables X\, . . . , X n assuming values in Gm- This becomes our state 
space, and the set of n random variables is a finite stochastic process. Its full characterization 



is in term of all finite dimensional distributions of the following kind (for 1 < m < n) 17] 



Px u ...,x m (xi, ...,x m )= P(Xi = xi, . . . ,X m = x m ), (1) 
where P(-) denotes the probability of an event with the values Xi running on all the possible 



graphs Qu of Gm- The finite dimensiona' 
compatibility conditions of Kolmogorov 



distributions defined in equation ([I]) obey the two 
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namely a symmetry condition 



Px u ...,x m (xi, ...,x m ) = Px^x^ix^, ■ ■ ■ ,x Vm ) (2) 

for any permutation (tti, . . . , ir m ) of the m random variables (this is a direct consequence of 
the symmetry property for the intersection of events) and a second condition 

Pxi,...,x m (x±, . . . , x m ) = ^ Px u ...,x m ,x m+1 {xi, - ■ ■ , X mi X m +i ) 

(3) 

in 

as a direct consequence of total probability. 

Among all possible stochastic processes on Gm, we will consider homogeneous Markov 
chains. They are fully characterized by the initial probability 

p(x)=p Xl (x)=¥(X 1 = x) (4) 

and by the transition probability 

P(x, y) = P(X m+1 = y\X m = x) (5) 

that does not depend on the specific value of m (hence the adjective homogeneous). Note 
that it is convenient to consider the initial probability as a row vector with 2 M( - M+1 ^ 2 entries 
with the property that 

E = !» ( 6 ) 

x£G M 

and the transition probability as a 2 M<yM+l ^ 2 x 2 M ( M+1 )/ 2 matrix, also called stochastic 
matrix with the property that 



For a homogeneous Markov chain, the finite dimensional distributions are given by 

PXi,...,X m (^l, • • • , X m ) = P(xi)P(xi, X 2 ) ■ ■ ■ P(x m -i, X m ). (8) 

It is a well known fact that the finite dimensional distributions in equation (JH]) do satisfy 
Kolmogorov's conditions ([2]) and ([3]). Kolmogorov's extension theorem then implies the 
existence of Markov chains 17|. Marginalization of equation (jSJ) leads to a formula for 



Px m {x m ) = ^{Xm = x m ), this is given by 

Pxjx m )= p(x 1 )P m -\x 1 ,x m ) (9) 

xieGM 

where P m ~ l (xi,x m ) is the entry (xi,x m ) of the (m — l)-th power of the stochastic matrix. 
Note that, from equation (Q and homogeneity one can prove the Markov semi-group property 

P m+r (x,y)= J2 P m (x,z)P r (x,y). (10) 

z&Gm 

Starting from the basic Markov process with the set of graphs as space state, we can also 
consider other auxiliary processes. Just to mention few among them, we recall: 

• the process counting the number of edges (i.e., the sum of the adjacency matrix A); 

• the process recording the degree of the graph (i.e., the marginal total of the adjacency 
matrix A); 

• the process which measures the cardinality of the strongly connected components of 
the graph. 

Notice that the function of a Markov chain is not a Markov chain in general, and therefore 
the study of such processes is not trivial. 

Under a more combinatorial approach, one can consider also the process recording the 
permanent of the adjacency matrix A. We recall that the permanent of the matrix A is 
given by 

M 

perm(A) = H^^) ( H ) 

where Sm is the symmetric group on the set {1, . . . , M}. The permanent differs from the 
best known determinant only in the signs of the permutations. In fact, 

M 

det(A)= X)(-l) w nV« ( 12 ) 



where |er| is the parity of the permutation a. Notice that the permanent is in general harder 
to compute than the determinant, as Gaussian elimination cannot be used. However, the 
permanent is more appropriate to study the structure of the graphs. It is known, see for 
instance [if]], that the permanent of the adjacency matrix counts the number of the bijective 
functions : Vm — > Vm- The bijective functions are known in this context as perfect 
matchings, i.e., the rearrangements of the vertices consistent with the edges of the graph. 
The relations between permanent and perfect matchings are especially studied in the case 
of bipartite graphs, see 18| for a review of some classical results. 

Moreover, we can approach the problem also from the point of view of symbolic compu- 
tation, and we introduce the permanent polynomial, defined for each adjacency matrix as 
follows. Let Y be an M x M matrix of variables Y = (2/jj)™ = i- The permanent polynomial 
is the polynomial 

pperm(A) = perm(Y A) , (13) 

where denotes the element-wise product. For example, the polynomial determinant of the 
adjacency matrix 
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1 1 
1 1 



introduced above is 



^ y 1>2 ?/i,3 ^ 



pperm(A) = det 



2/1,22/2,32/3,1 + 2/1,32/3,22/2,1 • 



2/2,1 2/2,3 
\^ 2/3,1 2/3,2 J 

The permanent polynomial in Equation (TlBl) is a homogeneous polynomial with degree M 
and it has as many terms as the permanent of A, all monomials are pure (i.e., with unitary 
coefficient) and each transition of the Markov chain from the adjacency matrix Ai to the 
matrix A 2 induces a polynomial pperm(A 2 ) — pperm(Ai). 

Finally, is is also interesting to consider conditional graphs. With this term we refer to 
processes on a subset of the whole family of graphs Gm- For instance we may require to move 
only between graphs with a fixed degree, i.e., between adjacency matrices with fixed row 
(and column) totals. In such a case, also the construction of a connected Markov chain in 
discrete time is an open problem, recently approached through algebraic and combinatorial 
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techniques based on the notion of Markov basis, see 19|-|2l). This research topic, named 
Algebraic Statistics for contingency tables, seems to be promising when applied to adjacency 
matrices of graphs. 



B. Ingredient 2: a semi-Markov counting process 

Let Ji, . . . , J n , . . . be a sequence of positive independent and identically distributed (i.i.d.) 
random variables interpreted as sojourn times between events in a point process. They are 
a renewal process. Let 

n 

T n = J2 J * ( 14 ) 

i=l 

be the epoch (instant) of the n-th event. Then, the process N(t) counting the events occurred 
up to time t is defined by 

N(t) = max{n : T n < t}. (15) 

A well-known (and well-studied) counting process is the Poisson process. If J ~ exp(A), 
one can prove that 

P(n,t) =F(N{t) =n)= exp(-At)^^. (16) 

The proof leading to the exponential distribution of sojourn times to the Poisson distribution 
of the counting process is rather straightforward. First of all one notices that the event 
{N(t) < n + 1} is given by the union of two disjoint events 

{N(t) < n + 1} = {N(t) <n}U {N(t) = n}, (17) 

therefore, one has 

P(JV(t) = n) = ¥(N(t) < n + 1) - P(JV(t) < n); (18) 

but, by definition, the event {N(t) < n} coincides with the event {T n > t}. Therefore, from 
equation (fl8l) . one derives that 

F(N(t) =n)= P(T n <t)- P(T„ +1 < t). (19) 



The thesis follows from equation (1141) . The cumulative distribution function of T n is the 
n-fold convolution of an exponential distribution, leading to the Erlang (or Gamma) distri- 
bution 



n-l .\fc 



P(T n <t) = l - ^exp(-At)^, (20) 

fc=0 
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and, by virtue of equation ( Fl9|) . the difference P(T n < t) — P(T n +i < t) gives the Poisson 
distribution of equation ffl6|) . Incidentally, it can be proved that N(t) has stationary and 
independent increments. 

One can also start from the Poisson process and then show that the sojourn times are 
i.i.d. random variables. The Poisson process can be defined as a non-negative-integer-valued 
stochastic process N(t) with N(0) = and with stationary and independent increments (i.e. 
a Levy process; well, it must also be stochastically continuous, that is it must be true that 
for all a > 0, and for all s > lim t ^ s P(|iV(£) — iV(s)| > a) = 0) such that its increment 
N(t) — N(s) with (0 < s < t) has the following distribution for n > 

F{N{t) - N{s) =n)= exp(-A(t - s)) ^ - S ^ . (21) 

Based on the definition of the process, it is possible to build any of its finite dimensional 
distributions using the increment distribution. For instance P(N(ti) = n 1; iV(t 2 ) = n 2 ) with 
t 2 > t\ is given by 

P(N(t 1 )=n u N(t 2 )=n 2 ) = 

= P(AT(t 1 ) = ni )P(N(h) - N(h) =m- m) 

= eM-^h)^ eM-Kh - tO) [Kt ; " tl)] Z ni ■ (22) 
nil (n 2 -nij! 

Every Levy process, including the Poisson process is Markovian and has the so-called strong 

Markov property roughly meaning that the Markov property is true not only for deterministic 

times, but also for random stopping times. Using this property, it is possible to prove 

that the sojourn times are independent and identically distributed. For N(0) = 0, let 

T n = inf{t : N(t) = n} be the n-th epoch of the Poisson process (the time at which the 

n-th jump takes place) and let = Tk — T^-i be the k-th sojourn time (T = 0). For what 

concerns the identical distribution of sojourn times, one has that 

P(Ti >t)= P(Jj >t)= P(N(t) = 0) = exp(-At), (23) 

and for a generic sojourn time Tk — Tk-i, one finds 

Tk - T k+ i = inf {t - T^ : N(t) = k} 

= inf {t - T^ : N(t) - iV(7W) = 1} 

= inf {t - T^! ■■ N(t - Tk-!) = 1, N(T k -!) = 0} 

= inf {t : N(t) = 1, N(0) = 0}, (24) 

9 



where = denotes equality in distribution and the equalities are direct consequences of the 
properties defining the Poisson process. The chain of equalities means that every sojourn 
time has the same distribution of J\ whose survival function is given in equation ( 123]) . As 
mentioned above, the independence of sojourn times is due to the strong Markov property. 
As a final remark, in this digression on the Poisson process, it is important to notice that 
one has that its renewal junction H(t) = f K(N{t)) is given by 

H{t) = Xt (25) 

i.e. the renewal function of the Poisson process linearly grows with time, whereas its renewal 
density h(t) defined as 

h(t) ^ ^ (26) 

is constant: 

hit) = A. (27) 

Here, for the sake of simplicity, we shall only consider renewal processes and the related 
counting processes (see equations (fl4l) and (f!5l) ). When sojourn times are non-exponentially 
distributed, the corresponding counting process N(t) is no longer Levy and Markovian, but 



it 



Delongs to the class of semi-Markov processes further characterized in the next section [2 



25|. If ipj(t) denotes the probability density function of sojourn times and Vj(t) = P(J > t) 



is the corresponding survival function, it is possible to prove the first renewal equation 

H(t) = l-Vj(t)+ f H(t-u)tljj(u)d Uj (28) 
Jo 



as well as the second renewal equation 

h(t) = ipj(t) + [ h(t-u)ipj(u)du. (29) 
Jo 

The second renewal equation is an immediate consequence of the first one, based on the 
definition of the renewal density h{t) and on the fact that tpj(t) = —d^j(t)/dt. The first 
renewal equation can be obtained from equation (fl9l) which is valid in general and not only 
for exponential waiting times. One has the following chain of equalities 

oo oo 

H(t) = E(N(t)) = J2^(N(t) =n) = J>(P(T n < t) - F(T n+1 < t)) 

n=0 n=0 

oo oo 

= ^P(T„<t) = ^F n (t), (30) 



n=l n=l 
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where F n (t) is the cumulative distribution function of the random variable T n , a sum of 
i.i.d. positive random variables. Let f n (t) represent the corresponding density function. By 
taking the Laplace transform of equation (I3"UI) and using the fact that 

f n (s) = $j(s)] n , (31) 

one eventually gets 



S — S — S z — ' S 1 — ?/) rfsl 

n=l n=l n=l m=0 1 rA 5 J 



or (as |V ; j('5)| < 1 for s 7^ 0) 



(32) 



[l-Ms))H(s) = ^-; (33) 



the inversion of equation ( 133|) yields the first renewal equation ( j28J) . 

If the sojourn times have a finite first moment (i.e. \xj = E(J) < 00), one has a strong 
law of large numbers for renewal processes 

lim = J_ a .s. (34) 

and as a consequence of this result, one can prove the so-called elementary renewal theorem 

lim = — . (35) 

The intuitive meaning of these theorems is as follows: if a renewal process is observed a 
long time after its inception, it is impossible to distinguish it from a Poisson process. As 
mentioned in section [J, the elementary renewal theorem can explain the ubiquity of the 
Poisson process. After a trasient period, most renewal processes behave as the Poisson 
process. However, there is a class of renewal processes for which the condition E( J) < 00 is 
not fulfilled. These processes never behave as the Poisson process. A prototypical example 
is given by the renewal process of Mittag-Leffler type introduced by one of us together with 



F. Mainardi and R. Gorenflo back in 2004 
will be given in one of the examples below. 
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27] . A detailed description of this process 
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C. Putting the ingredients together 



Let Xi, . . . , X n represent a (finite) Markov chain on the state space Gm, we now introduce 
the process Y(t) defined as follows 

Y(t)=X m , (36) 

that is the Markov chain X n is subordinated to a counting process N(t) coming from a 
renewal process as discussed in the previous subsection, with X n independent of N(t). In 
other words, Y(t) coincides with the Markov chain, but the number of transitions up to 
time t is a random variable ruled by the probability law of N(t) and the sojourn times in 
each state follow the law characterized by the probability density function if)j(t), or, more 
generally, by the survival function ^j(t). 
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As already discussed, such a process belongs to the class of semi-Markov processes 
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281 ]. i.e. for any A C Gm and t > we do have 

P(X n e A, J n < 1 1 X Q , . . . , X n _i, Ji, . . . , J n _i) = F(X n eA,J n <t\ X„_ x ) (37) 

and, if the state X n _i = x is fixed at time t n _i, the probability on the right-hand side will 
be independent of n. Indeed, by definition, given the independence between the Markov 
chain and the counting process, one can write 

F(X n e A, J n < 1 1 X , . . . , X n -!, Ji, . . . , Jn-i) = P(X n e A I X n _x = x)F{J n < t) 

= P(x,A)(l-*j(t)), (38) 

where 

P(x,A) = ^2P(x,y). (39) 
yeA 

Equation is a particular case of (13T1) . 

It is possible to introduce a slight complication and still preserve the semi-Markov prop- 
erty. One can imagine that the sojourn time in each state is a function of the state itself. In 
this case P(J„ < t) is no longer independent of the state of the random variable X n _i and 
equation (1351) is replaced by 

P(X n e A, J n < 1 1 X , . . . , X n _ 1; J u . . . , J n _i) = F(X n e A \ X n _, = x)F(J n < t\X n ^ = x) 

= P(x,A)(l-* x j(t)), (40) 

12 



where *&j(t) denotes the state-dependent survival function. However, in this case, the ran- 
dom variable T n is still the sum of independent random variables, but they are no-longer 
identically distributed, and the analysis of the previous section has to be modified in order 
to take this fact into proper account. 



III. EXAMPLES 



In order to show the behavior of the stochastic processes described in the previous sections 
we have simulated the distribution of two stopping times in two different situations. The 
simulations have been written in R, see 29] and the source files are available as supplementary 



online material. Notice that some specific packages for the analysis of graph structures are 



available, see for instance 



30j. However, we have used only the R-base commands, as our 



examples can be analyzed easily without any additional package. 

The examples in this section are designed to introduce the reader to the simulation 
algorithms in a framework as simple as possible. An extended example about a model of 
interbank market will be discussed in the next section. 

In our examples we use the Mittag-Leflier distribution for the sojourn times. We recall 
that the Mittag-Leflier distribution has survival function given by 



*j{t) = P( J >t) = E p (-t 



where Ep(z) is the one-parameter Mittag-Leflier function defined as 



(41) 



n=0 



r(n/3 + l)' 



(42) 



for < (3 < 1. There are two strong reasons for this choice. The first one is that many 



analytical resu 



son process 
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ts are available on the Mittag-Leflier renewal process a.k.a. fractional Pois- 
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3 lM33(| . The second reason is that the Mittag-Leflier distribution is the 



repeated-thinning limit of heavy-tailed sojourn-time distributions with algebraically decay- 
ing tails with exponent < < 1 [271 ] . For /3 = 1, the exponential distribution is recovered 
from (HE). 
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First example (Triangles) 
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M=10 nodes 

FIG. 1: Box-plot of the distribution of the stopping times with varying f3 for Example A. 
A. First example 

In this example we consider graphs without self-loops. Let us consider a fixed number 
M of vertices and define a process as follows: 

• At the time 0, there are no edges in the graph; 

• At each time, we choose an edge e with uniform distribution on the 2 — 2 ' edges. If 
e belongs to the graph we remove it; if e does not belong to the graph we add it; 

• The stopping time is defined as the first time for which a triangle appears in the graph. 

To simulate the distribution of the stopping times we have used 10, 000 replications. As 
the Mittag-Leffler distribution is heavy-tailed, the density plot and the empirical distribution 
function plot are not informative. Thus, we have reported the box-plot, to highlight the 
influence of the outliers. 

With a first experiment, we have studied the influence of the (3 parameter. In graph with 
M — 10 nodes, we have considered the sojourn times with a Mittag-Leffler distribution with 
different (3 parameter, namely j3 = 0.90.0.95, 0.98, 0.99. The box-plot are displayed in Figure 
[lj and some numerical indices are in Table [B 

Our results show that: 
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beta: 0.9 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

3.628e-01 8.010e+00 1.275e+01 3.146e+01 2.026e+01 5.475e+04 



beta: 0.95 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

3.249e-01 7.545e+00 1.144e+01 2.086e+01 1.689e+01 3.292e+04 



beta: 0.98 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

0.2607 7.2960 10.8600 12.9500 15.0500 2704.0000 



beta: 0.99 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

0.5373 7.2190 10.6300 12.3400 14.6700 2487.0000 

TABLE I: Summary statistics for Example A with varying f3. 

• the outliers are highly influenced from the value of f3. This holds, with a less strong 
evidence, also for the quartiles Ql and Q3; 

• the median is near constant, while the mean is affected by the presence of outliers. 

With a second experiment, we have considered a fixed parameter (3 = 0.99, but a variable 
number of vertices M ranging from 5 to 50 by 5. In Figure [2] we present the box-plots of 
the stopping time distribution and the trends of the mean and the median. 

From this graphs we can notice that: 

• the presence of outliers is more relevant in the graph with a large number of nodes; 

• the mean and the median are roughly linear, but the trend of the median looks more 
stable. 
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First example (Triangles) 



First example (Triangles) 
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FIG. 2: Box-plot (left), mean and median (right) of the distribution of the stopping times with 
varying M for Example A. 

B. Second example 

Let us consider a population with individuals {1, . . . , M} and suppose that the person 
labelled 1 has to share some information. At a first random time, he chooses another 
individual with random uniform probability and shares the information with him. At a 
second random time, one person who has the information chooses an individual among the 
other (M — 1) and shares again the information. Note that each individual shares the 
information with another one, no matters if he has already or not the information. In terms 
of graphs, we define a process as follows: 

• At the time 0, there are no edges in the graph; 

• At each time, we choose a vertex m connected with 1 and we add to the graph an edge 
among (m, 1), (m, 2), . . . , (m, m — 1), (m, m + 1), . . . , (m, M — 1), (m, M) with random 
uniform distribution. If the chosen edge is already in the graph we do nothing; 

• The stopping time is defined as the first time for which the whole graph is connected. 

The experimental settings for this example are the same as for Example A. With a fixed 
number of vertices M = 10 and varying {3 as above, we obtain the box-plots in Figure El 
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beta: 0.9 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

3.786e+00 2.154e+01 3.207e+01 8.811e+01 4.938e+01 2.715e+05 



beta: 0.95 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

3.393 19.410 27.050 41.550 38.260 12140.000 



beta: 0.98 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

3.565 18.230 24.980 33.530 33.970 19600.000 



beta: 0.99 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

4.738 17.690 23.940 27.160 32.310 1701.000 



TABLE II: Summary statistics for Example B with varying (3. 

and the numerical summary in Table ILT1 

From this results we can see that the outliers are highly influenced from the value of j3, 
while the variation of the quantiles Ql and Q3 is much lower. Also in this example, the 
mean is affected by the presence of outliers. 

With the second experiment with a variable number of vertices M ranging from 5 to 50 
by 5, we obtain the plots displayed in Figure HI The conclusions are the same as in the 
previous example. 



IV. EXTENDED EXAMPLE. AN INTERBANK MARKET 

In this section we present a simple model for interbank markets. It serves the purpose of 
illustrating the modelling potentialities of the ideas presented above. 

This example deals with an interbank market receiving loan requests from the corporate 
sector at random times. For the sake of readability, in this section we will use the symbol 
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Second example (Total connection) 



beta 
M=10 nodes 



FIG. 3: Box-plot of the distribution of the stopping times with varying (3 for Example B. 
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FIG. 4: Box-plot (left), mean and median (right) of the distribution of the stopping times with 
varying M for Example B. 



Atfc instead of for the k-th inter-arrival duration and we will denote the epochs at which 
loan requests are made with the symbol instead of In this more realistic example, we 
will briefly discuss the difficulties that must be faced when one desires to go beyond a mere 
phenomenological description of reality. 



18 



Assets 


Liabilities 


C\ n : liquidity 

L\ : loans to the corporate sector 
Ct : loans to other banks 


Df n : total (households' and firms') deposits 
-Bi : debt with other banks 

''Ti 

Et n : equity (net worth) 



TABLE III: Balance sheet entries of bank b at time t 



We consider an interbank market characterized by M banks that demand and supply 
liquidity at a given interest rate r#. Each bank b is described at any time by its balance 
sheet, as outlined in Table III II The market is decentralized and banks exchange liquidity 
by means of pairwise interactions. Banks lend money also to the corporate sector at the 
constant rate rc > re and all corporate and interbank loans are to be repayed after T 
units of time. We stipulate that the loan requests from the corporate sector to the banking 
system are the events triggering the interbank market and we model these events as a Poisson 
process of parameter A. In particular, we state that, at exponentially distributed intervals of 
time At n = t n — t n _i, a loan request of constant amount i is submitted from the corporate 
sector to a bank chosen at random with uniform distribution among the M banks. As in 
the previous examples, in principle, the Poisson process can be replaced by any suitable 
counting process. Let us denote the chosen bank with the index i and the time at which the 
loan is requested as t n . If C\ < £, the chosen bank is short of liquidity to grant the entire 
amount of the loan. Given the interest rate spread between rc and r#, the profit-seeking 
bank enters the interbank market in order to borrow at the rate rn the amount £ — CI 

L n — 1 

necessary to grant the entire loan. In the interbank market, a new bank is then chosen at 
random with uniform distribution among the remaining M — 1 banks. Let us denote with 
j the new chosen bank. If bank j has not enough liquidity to lend the requested amount, 
i.e., if Cj n _ 1 < I — C\ _ , then a new bank h is again chosen at random among the remaining 
M — 2 ones to provide the residual liquidity, and so on. This process in the interbank market 
continues till the liquidity amount I — C\ necessary to bank i is collected. 

Finally, as soon as the loan i is provided to the corporate sector, we stipulate that the 
deposits as well as the liquidity of any bank b, being b = 1, . . . , M, is increased by the amount 
u}^£, where u\ are random numbers constrained by ^2i b oj\ = 1. The rationale behind this 
choice is that a loan, when is taken and spent, creates a deposit in the bank account of the 
agent to whom the payment is made; for instance, when the corporate sector gets a loan to 
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pay wages to workers or to pay investments to capital goods producers, then the deposits at 
the M banks of the agents receiving the money are increased by a fraction of the borrowed 
amount I. We assume that these deposits are randomly distributed among the M banks. 

To give a clearer idea on how the balance sheets of banks evolve after an event in the 
interbank market, let us consider an example where at time t n the corporate sector requests 
a loan £ to the randomly selected bank i, which, being short of liquidity (i.e. C\ < £), 
needs to enter into interbank market where it borrows a loan of amount i — C\ from the 

t-n— 1 

randomly selected bank j. We suppose here > t — C\ , therefore no other lending 

banks enter the interbank market. According to the model outlined above, at the end of the 
interbank market session, the balance sheets of bank i and of bank j change as outlined in 
Table ED 
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TABLE IV: Dynamics of balance sheet entries of bank i (lender to the corporate sector and borrower 
in the interbank market) and bank j (lender in the interbank market) at time t n when both the 
corporate loan t and the related interbank loan £ — C\ n _ x are granted. 

Once the assets and the debt liabilities entries of any bank are updated following the 
lending activity to the corporate sector and the interbank market outcomes, the equity is 
then updated as residual according to the usual accounting equation: 

E h . = C h . + L h , + C h , - D h , - B\ . (43) 

It is worth noting that, as reported in Table IIVI the equity of both bank i and j does not 
change from t n -\ to t n . This result is obtained by computing the new equity levels at time 
t n using pB"]) but should not be a surprise given that lending and borrowing clearly change 
the balance sheet entries of banks but not their net worth at the time the loan is granted 
or received. Indeed, the net worth of the lending banks is increased by the interest revenues 
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when the corporate loan as well as the interbank loan is repaid together with the interest 
amounts. In particular, equity of bank i is increased by rci — tb (i — CI _ ), while equity 
of bank j is increased by {£ — C\ ). Table |V] shows how balance sheet entries change 
at time t m = t n + T when the two loans are paid back. It is worth noting again that the 
equity dynamics is consistent with the dynamics of other balance sheet entries, according 
to ( 143|) . Finally, as granting a bank loan to the corporate sector increases private deposits 
at banks, also the opposite holds when a loan is paid back. The repayment of the loan I 
together with interests rc t corresponds to a reduction of private deposits, as well as of the 
related banks' liquidity, of the same amount. As in the previous case, we assume that the 
reduction (1 +rc) £ is uniformly and randomly distributed among the M banks with weights 
bj h , , where b = 1, . . . , M. 
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TABLE V: Dynamics of balance sheet entries of bank i (lender to the corporate sector and borrower 
in the interbank market) and bank j (lender in the interbank market) at time t m = t n + T when 
both the corporate loan I and the related interbank loan I — C\ are paid back. 

We can then define a M x M adjacency matrix A representing the graph associated to 
the interbank market, where the nodes of the graph correspond the M banks and the edges 
to the lending and borrowing relationships among banks. At variance with the previous 
discussion and examples, here, it is meaningful to consider directed graphs and therefore 
the matrix can be asymmetric. In particular, if bank j is lending money to bank i, we set 
Aj t i = 1, but we may have Aij = 1 or Aij = 0, depending if bank i is lending or not money 
to bank j. The situation where both Ajj and Aij are set to 1 is not contradictory but it 
means that two loans have been granted in the two opposite directions, i.e. from bank i to 
bank j and from bank j to bank i, at different times. The time evolution of the adjacency 
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matrix depends on the evolution of events in the interbank market. In particular, when the 
first loan from bank j to bank i is paid back, Ajj is again set to 0, provided that no new 
loans have been granted by bank j to bank i in the meantime, if this happens the value of 
Aj t i remains at 1 till there are debts of bank i to bank j. If this is required by the particular 
application, it is even possible to consider weighted graphs where the entry A it j contains the 
value of the loan from bank % to bank j. 

The dynamics in the interbank market can then be represented as a Markov chain on 
graphs subordinated to the Poisson process representing the random events of loan requests 
to the banking system by the corporate sector. It is worth noting that the Markov process 
and the Poisson process are independent here, however, the transition probabilities of the 
Markov process are not fixed ex ante but depends on the endogenous evolution of the balance 
sheets of banks. Therefore, here, the Markov process is not homogeneous. 



V. CONCLUDING CONSIDERATIONS 



We have discussed a model of graph dynamics based on two ingredients. The first in- 
gredient is a Markov chain on the space of possible graphs. The second ingredient is a 
semi-Markov counting process of renewal type. The model consists in subordinating the 
Markov chain to the semi-Markov counting process. In simple words, this means that the 
chain transitions occur at random time instants called epochs. This model takes into account 
the fact that social interactions are intrinsically volatile and not permanent. 

Note that state dependent subordination (see equation (J1D]) ) gives rise to very interesting 
dynamics from the ergodicity viewpoint 34(. In order to illustrate this fact, let us consider 
a simple two-state aperiodic and irreducible Markov chain with the following transition 
probability matrix: 

' 0.1 0.9 
0.9 0.1 

In this case, the invariant measure is uniform and it is given by 



P 



P 



1 1 

2' 2 



meaning that the probability of finding each state at equilibrium is 1/2. Now, let us call A 
the first state and B the second state. Let the sojourn time in A be exponentially distributed 
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with parameter Xa and the sojourn time in B still exponentially distributed with parameter 
Xb- If a single realization of this process is considered, the average time of permanence in 
state A will be given by 1/Xa and the average time of permanence in B will be given by 
1/Xb- Therefore, if Xa^ Xb, then the ratio of average sojourn times will be different from 
1. In other words, for this simple model, the fraction of sojourn times is not equal to the 
fraction of the ensemble measure: a signal of non-ergodicity. 

Finally, with reference to the examples discussed above, this kind of modeling can be 
used for risk evaluation. Given a loss function, a function that gives the losses when adverse 
events take place, the risk function is defined as the expected value of the loss function. 
With our approach, one can derive the probability of the adverse events as a function of 
time and use this measure to evaluate the risk function. 

VI. ACKNOWLEDGEMENTS 

FR wishes to acknowledge an INDAM (Istituto Nazionale di Alta Matematica, Italy) 
grant with which this work was partially funded. Long ago (in 2005 and 2006), ES discussed 
these topics with Gilles Daniel, Lev Muchnik and Sorin Solomon. He now gratefully ac- 
knowledges these discussions. ES wishes to thank Elena Akhmatskaya for recent discussion 
on these issues. MR and ES first discussed these issues during a visit of ES to Reykjavik 
University sponsored by the ERASMUS program. ES is also grateful to Universitat Jaume 
I for the financial support received from their Research Promotion Plan 2010 during his 
scientific visit in Castellon de la Plana where parts this paper were written. 



[1] D.J. Watts and S.H. Strogatz, Collective dynamics of 'small world' networks, Nature 393, 
440-442, 1998. 

[2] R. Durrett, Random graph dynamics, Cambridge University Press, Cambridge UK, 2007. 
[3] A.L. Barabasi and R. Albert, Emergence of scaling in random networks, Science 286, 509-512, 
1999. 

[4] S. Milgram, The small world problem, Physiology Today 2, 60-67, 1967. 
[5] M. Granovetter, The strength of weak ties, American Journal of Sociology 78, 1360-1380, 
1973. 

23 



[6] G.U. Yule, A mathematical theory of evolution, based on the conclusions of Dr. J.C. Willis, 
Philosophical Transactions of the Royal Society of London B 213, 21-87, 1925. 

[7] H.A. Simon, On a class of skew distribution functions, Biometrika 42, 425-440, 1955. 

[8] U. Garibaldi and E. Scalas, Finitary probabilistic methods in econophysics, Cambridge Uni- 
versity Press, Cambridge UK, 2010. 

[9] A.L. Barabasi, Linked. How everything is connected to everything else and what it means for 
business, science, and everyday life, Plume, New York NY, 2003. 
[10] http: //en. wikipedia. org/wiki/Small_world_experiment 

[11] R.A. Howard, Dynamic probabilistic systems. Volume II: semi-Markov and decision processes, 

John Wiley & Sons, New York NY, 1971. 
[12] E. Scalas, R. Gorenflo, H. Luckock, F. Mainardi, M. Mantelli, and M. Raberto, Anomalous 

waiting times in high-frequency financial data, Quantitative Finance 4, 695-702, 2004. 
[13] E. Scalas, T. Kaizoji, M. Kirchler, J. Huber and A. Tedeschi, Waiting times between orders 

and trades in double- auction markets, Physica A 366, 463-471, 2006. 
[14] A.-L. Barabasi, Bursts: The Hidden Pattern Behind Everything We Do, Dutton Penguin, 

Boston, 2010. 

[15] F.R.K. Chung, Spectral Graph Theory, CBMS Regional Conference Series in Mathematics, 
1997. 

[16] N. Biggs, Algebraic Graph Theory, Second edition, Cambridge University Press, Cambridge, 
1993. 

[17] P. Billingsley, Probability and Measure, Wiley, New York, 1986. 

[18] A. Schrijver, Matching, Edge-Colouring, and Dimers, Graph-Theoretic Concepts in Computer 

Science, Lecture Notes in Computer Science, 2880, 13-22, 2003. 
[19] M. Drton, B. Sturmfels and S. Sullivant, Lectures on Algebraic Statistics, Birkhauser, Basel, 

2009. 

[20] F. Rapallo, Markov Bases and Structural Zeros, J. Symbolic Comput. 41, 164-172, 2006. 
[21] F. Rapallo and R. Yoshida, Markov bases and subbases for bounded contingency tables, Ann. 

Inst. Statist. Math. 62, 785-805, 2010. 
[22] E. Qinlar, Introduction to Stochastic Processes, Prentice-Hall, Englewood Cliffs, 1975. 
[23] O. Flomenbom, J. Klafter, Closed-Form Solutions for Continuous Time Random Walks on 

Finite Chains, Phys. Rev. Lett. 95, 098105, 2005. 



24 



[24] O. Flomenbom, R. J. Silbey, Path-probability density functions for semi-Markovian random 

walks, Phys. Rev. E 76, 041101, 2007. 
[25] J. Janssen and R. Manca, Semi-Markov Risk Models for Finance, Insurance and Reliability, 

Springer, New York, 2007. 
[26] E. Scalas, R. Gorenflo and F. Mainardi, Uncoupled continuous-time random walks: Solution 

and limiting behavior of the master equation, Phys. Rev. E 69, 011107, 2004. 
[27] F. Mainardi, R. Gorenflo and E. Scalas, A fractional generalization of the Poisson process, 

Vietnam Journal of Mathematics, 32 (SI), 53-64, 2004. 
[28] G. Germano, M. Politi, E. Scalas and R.L. Schilling, Stochastic calculus for uncoupled 

continuous-time random walks, Phys. Rev. E 79, 066102, 2009. 
[29] R Development Core Team, R: A Language and Environment for Statistical Computing, R 

Foundation for Statistical Computing, Vienna, Austria, 2010. 
[30] R. Gentleman, Elizabeth Whalen, W. Huber and S. Falcon, The 'graph' package, available at 

http://cran.r-project.org/, version 1.26.0, 2010. 
[31] D. Fulger, E. Scalas, and G. Germano, Monte Carlo simulation of uncoupled continuous-time 

random walks yielding a stochastic solution of the space-time fractional diffusion equation, 

Physical Review E 77, 021122, 2008. 
[32] L. Beghin and E. Orsingher, Fractional Poisson processes and related planar random motions, 

Electronic Journal of Probability 14, 1790-1826, 2009. 
[33] M.M. Meerschaert, E. Nane, and P. Vellaisamy, The fractional Poisson process and the inverse 



stable subordinator, http://arxiv.org/abs/1007.5051, 2010. 
[34] A. Saa and R. Venegeroles, Ergodic transitions in continuous-time random walks, Phys. Rev. 
E 82, 031110, 2010. 



25 



